#include "stdafx.h"
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <string.h> // memcpy
#include <fstream>
#include <omp.h>
#include "portable_time.h"

#ifndef M_PI
#define M_PI       3.14159265358979323846
#endif

double t_begin, t_end, time_counter;

double sigma = 1.0;
double k = 1.0;

double TStart, XStart, YStart, ZStart;
double TEnd, XEnd, YEnd, ZEnd;

double** UX;
double** UY;
double** UZ;

//ВРЕМЯ
double GetTime()
{
	return time_counter;
}

double f_begin(double x, double y, double z){ // начальное условие, функция u(0,x,y,z)=ksi(x,y,z)
	double f=0.0;
	if ((x<5.0) && (y<5.0) && (z<5.0))
	{
		f=sin(x*M_PI/5.0)+sin(y*M_PI/5.0)+sin(z*M_PI/5.0);
	}
	return f;
}

double f(double t, double x, double y, double z){ // функция, задающая внешнее воздействие
	double fxyz=0.0;
	if (((x>=5.0) && (x<=5.2)) && ((y>=5.0) && (y<=5.2)) && ((z>=5.0) && (z<=5.2)))
		fxyz=4.0;
	return fxyz;
}

double f_left_x(double time, double y, double z){ // левое граничное условие по x
	return 0;
}

double f_right_x(double time, double y, double z){ // правое граничное условие по x
	return 0;
}

double f_left_y(double time, double x, double z){ // левое граничное условие по y
	return 0;
}

double f_right_y(double time, double x, double z){ // правое граничное условие по y
	return 0;
}
double f_left_z(double time, double x, double y){ // левое граничное условие по z
	return 0;
}

double f_right_z(double time, double x, double y){ // правое граничное условие по z
	return 0;
}

void calculation(	
	int Nx, int Ny, int Nz,	//количество интервалов
	int NTimes     //количество интервалов времени
	)
{
	t_begin = PortableGetTime();
	//
	double XStart = 0.0; double YStart = 0.0; double ZStart = 0.0;
	double TStart = 0.0;
	/*
	double deltaX = (XEnd - XStart) / Nx;
	double deltaY = (YEnd - YStart) / Ny;
	double deltaZ = (ZEnd - ZStart) / Nz;
	double deltaT = (TEnd - TStart) / NTimes;
	*/

	//сколько точек
	double ***data[2];
	int stepsX = Nx + 1;
	int stepsY = Ny + 1;
	int stepsZ = Nz + 1;
	int stepsT = NTimes + 1;

	//пусть интервалов столько же, сколько точек
	double deltaX = (XEnd - XStart) / stepsX;
	double deltaY = (YEnd - YStart) / stepsY;
	double deltaZ = (ZEnd - ZStart) / stepsZ;
	double deltaT = (TEnd - TStart) / stepsT;

	data[0] = new double**[stepsX];
	data[1] = new double**[stepsX];
	for (int i=0; i<stepsX; i++)
	{
		data[0][i] = new double*[stepsY];
		data[1][i] = new double*[stepsY];
		for (int j=0; j<stepsY; j++)
		{
			data[0][i][j] = new double[stepsZ]();
			data[1][i][j] = new double[stepsZ]();
		}
	}

	UX = new double*[stepsT];
	UY = new double*[stepsT];
	UZ = new double*[stepsT];
	for (int i=0; i<stepsT; i++)
	{
		UX[i] = new double[stepsX];
		UY[i] = new double[stepsY];
		UZ[i] = new double[stepsZ];
	}

	int projX, projY, projZ; // моменты времени, по которым мы делаем сечения
	projX = (stepsX - 1) / 2;
	projY = (stepsY - 1) / 2;
	projZ =(stepsZ - 1) / 2;

	// считаем в начальное время:
	for (int i=1; i<stepsX-1; i++)
		for (int j=1; j<stepsY-1; j++)
			for (int k=1; k<stepsZ-1; k++)
				data[0][i][j][k] = f_begin(XStart + i*deltaX, YStart + j*deltaY, ZStart + k*deltaZ);

	//пишем результаты для t=0
	for (int i=0; i<stepsX; i++)
		UX[0][i] = data[0][i][projY][projZ];
	for (int j=0; j<stepsY; j++)
		UY[0][j] = data[0][projX][j][projZ];
	for (int k=0; k<stepsZ; k++)
		UZ[0][k] = data[0][projX][projY][k];

	//считаем во все остальные моменты времени
	//for(int t = 1; t < stepsT; t++)
	//пусть моментов времени столько же, сколько и интервалов
	for(int t = 1; t < NTimes; t++)
	{
		// строим электронные плотности во все моменты времени
		int curr = t % 2;	// вычисление в n+1 момент времени
		int prev = (t-1) % 2; // по n-ному моменту времени

		int i;
#pragma omp parallel for shared(data) private(i)

		for (i=1; i<stepsX-1; i++)
		{
			for (int j=1; j<stepsY-1; j++)
			{
				for (int k=1; k<stepsZ-1; k++)
				{
					
					data[curr][i][j][k] = deltaT * (
						(data[prev][i+1][j][k] - 2*data[prev][i][j][k] + data[prev][i-1][j][k]) / (deltaX*deltaX)
						+
						(data[prev][i][j+1][k] - 2*data[prev][i][j][k] + data[prev][i][j-1][k]) / (deltaY*deltaY)
						+
						(data[prev][i][j][k+1] - 2*data[prev][i][j][k] + data[prev][i][j][k-1]) / (deltaZ*deltaZ)
						+
						f(TStart+(t-1)*deltaT, XStart+i*deltaX, YStart+j*deltaY, ZStart+k*deltaZ)
						-
						data[prev][i][j][k])
						+
						data[prev][i][j][k];

					/*
					data[curr][i][j][k] = (data[prev][i+1][j][k] - 2*data[prev][i][j][k] + data[prev][i-1][j][k]) / (deltaX*deltaX);
					data[curr][i][j][k] += (data[prev][i][j+1][k] - 2*data[prev][i][j][k] + data[prev][i][j-1][k]) / (deltaY*deltaY);
					data[curr][i][j][k] += (data[prev][i][j][k+1] - 2*data[prev][i][j][k] + data[prev][i][j][k-1]) / (deltaZ*deltaZ);

					data[curr][i][j][k] *= sigma;
					data[curr][i][j][k] -= k*data[prev][i][j][k];
					data[curr][i][j][k] += f(TStart+(t-1)*deltaT, XStart+i*deltaX, YStart+j*deltaY, ZStart+k*deltaZ);
					data[curr][i][j][k] *= deltaT;

					data[curr][i][j][k] += data[prev][i][j][k];
					*/
				};
			};
		};
		/*
		if (t==33)
		{
		FILE* testnums = fopen("file.txt", "w");
		for (i=1; i<stepsX-1; i++)
		for (int j=1; j<stepsY-1; j++)
		for (int k=1; k<stepsZ-1; k++)
		fprintf(testnums, "%lf(%i,%i,%i)\n", data[curr][i][j][k],i,j,k);
		fclose(testnums);
		}
		*/
		//собираем данные
		for (int i=0; i<stepsX; i++)
			UX[t][i] = data[curr][i][projY][projZ];
		for (int j=0; j<stepsY; j++)
			UY[t][j] = data[curr][projX][j][projZ];
		for (int k=0; k<stepsZ; k++)
			UZ[t][k] = data[curr][projX][projY][k];

	}

	for (int i=0; i<stepsX; i++)
	{
		for (int j=0; j<stepsY; j++)
		{
			delete [](data[0][i][j]);
			delete [](data[1][i][j]);
		};

		delete [](data[0][i]);
		delete [](data[1][i]);
	};
	//
	t_end = PortableGetTime();
	time_counter = t_end - t_begin;
	delete [](data[0]);
	delete [](data[1]);
}

int run(int argc, char **argv)
{
	XStart = YStart = ZStart = TStart = 0.0;

	XEnd = atof(argv[1]); // конец отрезка по x
	YEnd = atof(argv[2]); // конец отрезка по y
	ZEnd = atof(argv[3]); // конец отрезка по z

	int Nx = atoi(argv[4]); // количество интервалов, например 100
	int Ny = atoi(argv[5]); // количество интервалов, например 100
	int Nz = atoi(argv[6]); // количество интервалов, например 100

	double deltaX = (XEnd - XStart) / Nx;
	double deltaY = (YEnd - YStart) / Ny; 
	double deltaZ = (ZEnd - ZStart) / Nz; 

	TStart = 0.0;
	TEnd = atof(argv[7]); // конец времени моделирования
	int NTimes = atoi(argv[8]); // количество интервалов по времени, например 100

	double deltaT = (TEnd - TStart) / NTimes;

	calculation(Nx, Ny, Nz, NTimes);

	//вывод в файл
	FILE* filex=fopen("filex.txt","w");
	FILE* filey=fopen("filey.txt","w");
	FILE* filez=fopen("filez.txt","w");

	fprintf(filex,"%e\n", deltaX);
	fprintf(filey,"%e\n", deltaY);
	fprintf(filez,"%e\n", deltaZ);

	//
	for (int t = 0; t < NTimes;  t++)
	{
		fprintf(filex, "time = %i\n", t);
		fprintf(filey, "time = %i\n", t);
		fprintf(filez, "time = %i\n", t);
		for (int i = 0; i <= Nx; i++)
			fprintf(filex, "%lf ", UX[t][i]);
		for (int j = 0; j <= Ny; j++)
			fprintf(filey, "%lf ", UY[t][j]);
		for (int k = 0; k <= Nz; k++)
			fprintf(filez, "%lf ", UZ[t][k]);

		fprintf(filex,"\n");
		fprintf(filey,"\n");
		fprintf(filez,"\n");
	}

	fclose(filex);
	fclose(filey);
	fclose(filez);



	return 0;
};

double** getUX()
{
	return UX;
}

double** getUY()
{
	return UY;
}

double** getUZ()
{
	return UZ;
}

int  getMaxTime()
{
	return 100;
}

int getMaxCoord()
{
	return 100;
}